

#Getting the local path to R libraries
if(Sys.info()[["sysname"]]=="Windows"){
	library(Biobase)
}else{
	locFile =file("../Misc/localization_files.txt","r")
	for(line in readLines(locFile)){
		if(substr(line,0,9) == "Rlib_path"){
			Rlib_path = gsub(" ","",sub("=","",sub("Rlib_path","",line)))
		}
	}
	library(Biobase, lib.loc=Rlib_path)
}



expressionFile="../TXT/explvl.txt"
clinical_dataFile="../TXT/clinical_var.txt"
annotationFile="../TXT/annotation.txt"
tissue_nameFile="../TXT/tissuetype.txt"



#REFORMING FROM PHP
#loading data
expression<-read.table(expressionFile,sep="\t",row.names=1,header=TRUE)
colnames(expression)<-sub("^X","",colnames(expression))

clinical_data<-read.table(clinical_dataFile,sep="\t",row.names=1,header=TRUE)
colnames(clinical_data)<-sub("^X","",colnames(clinical_data)) 

annotation<-read.table(annotationFile,sep="\t",row.names=1,header=TRUE)
tissue_name<-read.table(tissue_nameFile,stringsAsFactors=FALSE)[,1]

#forming ExpressionSet
expressionset<- new("ExpressionSet",exprs = t(as.matrix(expression)),experimentData = new("MIAME",title=tissue_name))
pData(expressionset)<-clinical_data
if(!all(sort(featureNames(expressionset)) == sort(rownames(annotation))))stop("annotation rownames and expression rownames were not exactly equal")


#
#correlations<-data.frame("p-value"=vector())
#colnames(correlations)<-"p-value"
#for(clinVar in colnames(pData(expressionset))){
#	if(!class(expressionset[[clinVar ]])%in%c("factor","integer","numeric"))stop("clinVar must be either factor, integer or numeric when read in")
#	
#	if(class(expressionset[[clinVar ]])%in% "numeric"){
#		expressionset[[clinVar ]] <- as.numeric(expressionset[[clinVar ]])
#	}
#	clinData <- expressionset[[clinVar]]
#	for(probeset in featureNames(expressionset)){
#		expression <- exprs(expressionset)[probeset,]
#		model <- lm(expression ~ clinData)
#		if(class(expressionset[[clinVar ]])%in%"factor"){
#			correlations[paste(clinVar ,probeset),"p-value"]<-anova(model)["clinData","Pr(>F)"]
#		}
#		if(class(expressionset[[clinVar ]])%in%"numeric"){
#			correlations[paste(clinVar ,probeset),"p-value"] <- summary(model)[["coefficients"]]["clinData","Pr(>|t|)"]
#		}
#	}
#}
#correlations<-correlations[order(correlations[,"p-value"]),]
#
#









probesets<-featureNames(expressionset)

#other variables that don't need changing
horizontalScatterSpacing=0.02
doStatistics=TRUE
verbose=FALSE





plot_gene_vs_clinical<-function(expressionset,probesets,label,annotation=NULL,horizontalScatterSpacing=0.02,doStatistics=TRUE,verbose=TRUE,cexFactoring=1){
	for(label in colnames(pData(expressionset))){
		#ORIGINAL SCRIPT STARTS HERE
		#original arguments were: (expressionset,probesets,label,annotation=NULL,horizontalScatterSpacing=0.02,doStatistics=TRUE,verbose=TRUE){
		#loading a support function first:
		fun_plot_groupwise_expression_data_20100830<-function(groups,main="",ylab="expression",xlim=NULL,ylim=NULL,pointCol="black",pCutoff=0.05,vLineMultiplier=0.03,cexPValues=1,plotAxis=TRUE,plotN=TRUE,combinations=combn(names(groups),2),type="boxplot",horizontalScatterSpacing=0.1,groupXPosition=seq(1,length(groups)),cex.axis=1){
			#function to analyse different groups of data. Takes a named list of values "groups" and calculated
			#main						character with the plot title
			#ylab						the label of the y-axis
			#xlim						two numeric vector specifiying the limits of the x-axis. If left as NULL it tries to guess the best.
			#ylim						two numeric vector specifiying the limits of the y-axis. If left as NULL it tries to guess the best.
			#vLineMultiplier			is a numeric that can be used to control the vertical spacing between p-value lines
			#cexPValues					is a numeric for the size of the text in pValues
			#plotAxis					a logical indicating if the axis should be plotted
			#plotN						a logical indicating if the number of entries in each groups should be plotted with the axis
			#combinations				combinations to test with student's t-test. If NULL, this is omitted
			#type						character either "boxplot" or "dotplot"
			#horizontalScatterSpacing	a number - only used for dotplots. Decide how far the dots will be scattered horizontally
			#groupXPosition				a vector indicating the x-axis placement of groups. Can be used to specify custom breaks or spacing.
			#pointCol					A character vector specifying a colour for the points - defaults to black
			
			if(class(type)!="character")stop(paste("type must be of class character, not",class(type)))
			if(length(type)!=1)stop("type must be of length 1")
			if(!type%in%c("boxplot","dotplot"))stop("type must be either: 'boxplot' or 'dotplot'")
			if(class(groups)!="list")stop(paste("groups must be of class list, not",class(groups)))
			for(entry in names(groups)){
				if(!class(groups[[entry]])%in%c("NULL","numeric"))stop(paste("All list entries must be numeric vectors and",entry,"was not"))
			}
			if(type=="dotplot"){
				if(class(horizontalScatterSpacing)!="numeric")stop(paste("horizontalScatterSpacing must be of class numeric, not",class(horizontalScatterSpacing)))
				if(length(horizontalScatterSpacing)!=1)stop("horizontalScatterSpacing must be of length 1")
			}
			if(!class(groupXPosition)%in%c("integer","numeric"))stop(paste("groupXPosition must be of class numeric or integer, not",class(groupXPosition)))
			if(length(groupXPosition)!=length(groups))stop("groupXPosition must be of the same length as the number of groups")
			
			if(!pointCol%in%colors())stop(paste("Didn't recognize the colour",pointCol))
			
			vLine <- (max(unlist(groups),na.rm=TRUE)-min(unlist(groups),na.rm=TRUE)) * vLineMultiplier
			
			
			if(is.null(ylim)){
				
				ylim<-c(min(unlist(groups),na.rm=TRUE),max(unlist(groups),na.rm=TRUE)+10*vLine)
			}else{
				if(class(ylim)!="numeric")stop(paste("ylim must be of class numeric, not",class(ylim)))
				if(length(ylim)!=2)stop("ylim must be of length 2")
			}
			if(is.null(xlim)){
				xlim<-c(1-(length(groups)*0.05),length(groups)*1.05)
			}else{
				if(class(xlim)!="numeric")stop(paste("xlim must be of class numeric, not",class(xlim)))
				if(length(xlim)!=2)stop("xlim must be of length 2")
			}
			
			if(type=="boxplot"){
				categoriesCharacter<-vector()
				expression<-vector()
				for(level in names(groups)){
					categoriesCharacter<-c(categoriesCharacter,rep(level,length(groups[[level]])))
					expression<-c(expression,groups[[level]])
				}
				categories<-factor(x=categoriesCharacter,levels=names(groups))
				plot.default(
						NULL,
						xaxt="n",
						main=main,
						ylab=ylab,
						ylim=ylim,
						xlim=xlim,
						xlab=""
				)
				boxplot(expression ~ categories,
						add=TRUE,
						xaxt="n",
						at=groupXPosition,
						xlab="")
			}
			if(type=="dotplot"){
				
				plot(
						NULL,
						xaxt="n",
						ylab=list(ylab,cex=cexFactoring),
						xlab="",
						xlim=xlim,
						ylim=ylim,
						main=main
				)
				NumberOfDotsOfCexOneInOneLane<-50
				for(x in 1:length(groups)){
					level<-names(groups)[x]
					if(length(groups[[level]])>0){
						distribution<-hist(groups[[level]],plot=FALSE,breaks=NumberOfDotsOfCexOneInOneLane)
						bins<-c(ylim[1],distribution[["breaks"]],ylim[2])
						for(i in 1:(length(bins)-1)){
							ys<-groups[[level]][bins[i] < groups[[level]] & groups[[level]] < bins[i+1]]
							if(length(ys)>0){
								for(j in 1:length(ys)){
									xHere<-groupXPosition[x]-(length(ys)/2)*horizontalScatterSpacing + j*horizontalScatterSpacing - horizontalScatterSpacing/2
									points(y=ys[j],x=xHere,col=pointCol)
								}
							}
						}
					}
				}
			}
			
			
			#making the axis
			if(plotAxis){
				axisLabels<-vector()
				for(level in names(groups)){
					if(plotN){
						axisLabels<-c(axisLabels,paste("\n",level,"\nn=",length(groups[[level]]),sep=""))
						axisPadj=0.5
					}else{
						axisLabels<-c(axisLabels,level)
						axisPadj=0
					}
				}
				axis(1,at=groupXPosition,labels=axisLabels,padj=axisPadj,cex.axis=cex.axis*cexFactoring)
			}
			
			if(!is.null(combinations)){
				if(class(combinations)!="matrix")stop(paste("combinations must be of class matrix, not",class(combinations)))
				if(nrow(combinations)!=2)stop("combinations must have two lines")
				if(!all(unique(c(combinations[1,],combinations[2,]))%in%names(groups)))stop("One or more of the entries in combinations were not found in the names of groups")
				significantCount<-0
				for(i in 1:ncol(combinations)){
					p<-try(signif(t.test(groups[[combinations[1,i]]],groups[[combinations[2,i]]])[[3]],3),silent=TRUE)
					if(class(p)!="try-error"){
						if(p<pCutoff){
							xPos<-c(groupXPosition[match(combinations[1,i],names(groups))],groupXPosition[match(combinations[2,i],names(groups))])
							yMax<-ylim[2]-vLine-significantCount * vLine*3
							lines(x=xPos,y=c(yMax,yMax))
							lines(x=c(xPos[1],xPos[1]),y=c(yMax,yMax-vLine ))
							lines(x=c(xPos[2],xPos[2]),y=c(yMax,yMax-vLine ))
							text(x=mean(c(xPos[1],xPos[2])),y=yMax+vLine,labels=paste("P =",p),cex=cexPValues*cexFactoring)
							significantCount<-significantCount+1
						}
					}
				}
			}
		}
		
		
		
		if(class(expressionset)[1]!="ExpressionSet")stop(paste("expressionset must be of class ExpressionSet and not",class(expressionset)))
		if(class(probesets)!="character")stop(paste("probesets must be of class character and not",class(probesets)))
		if(class(label)!="character")stop(paste("label must be of class character and not",class(label)))
		if(length(label)!=1)stop(paste("label must be of length 1"))
		if(!label%in%colnames(pData(expressionset)))stop("The given label was not found in the expressionset")
		if(!class(expressionset[[label]])%in%c("factor","numeric","integer","Date"))stop(paste("The label",label,"was of class",class(expressionset[[label]]),"- it must be either factor, numeric, integer or Date"))
		if(!is.null(annotation)){
			if(!all(c("genename","genesymbol")%in%colnames(annotation))){
				stop("The given annotation file missed either the genename column or the genesymbol column")
			}
			
		}else{
			annotation<-data.frame(row.names=vector(),genesymbol=vector(),genename=vector())
		}
		
		if(verbose){print(paste("Plotting: now plotting",length(probesets),"from the set",experimentData(expressionset)@title,"at",label))}
		
		for (probeset in probesets){
			#FACTORS
			if(class(expressionset[[label]])%in%"factor"){
				groups<-list()
				for(group in levels(pData(expressionset)[,label])){
					groups[[group]]<-exprs(expressionset)[probeset,pData(expressionset)[,label]%in%group]
					
				}
				if(doStatistics){
					combinations=combn(names(groups),2)
				}else{
					combinations=NULL
				}
				
				if(length(levels(expressionset[[label]]))>6){
					cex.axis=0.4
				}else{
					cex.axis=0.7
				}
				
				
				fun_plot_groupwise_expression_data_20100830(
						groups=groups,
						main="",
						ylab="expression",
						xlim=NULL,
						ylim=NULL,
						pCutoff=0.05,
						vLineMultiplier=0.03,
						cexPValues=0.7,
						plotAxis=TRUE,
						plotN=FALSE,
						combinations=combinations,
						type="dotplot",
						cex.axis=0.7,
						horizontalScatterSpacing=horizontalScatterSpacing)
			}
			
			#NUMBERS and DATES
			if(class(expressionset[[label]])%in%c("integer","numeric","Date")){
				
				plot.default(
						type="p",
						y=exprs(expressionset)[probeset,],
						ylim=NULL,
						x=pData(expressionset)[,label],
						ylab="expression",
						xlab="",
						main="",
						cex=1
				)
				if(doStatistics & !class(expressionset[[label]])%in%"Date"){
					correlation<-signif(cor.test(exprs(expressionset)[probeset,],expressionset[[label]])[[4]],2)
					mtext(paste("Pearson correlation coefficient:",correlation),padj=7,side=1,cex=0.7*cexFactoring,adj=0)
				}
			}
			
			mtext(paste("Category:",label),padj=-1.4,cex=0.7*cexFactoring,adj=0)
			mtext(paste("Data set: ", expressionset@experimentData@title,sep=""),padj=-0.2,cex=0.7*cexFactoring,adj=0)
			
			
			if(probeset%in%rownames(annotation)){
				genedescription<-paste(probeset,annotation[probeset,"genesymbol"],"-",annotation[probeset,"genename"])
			}else{
				genedescription<-probeset
			}
			
			mtext(genedescription,side=1,cex=0.7*cexFactoring,adj=0,padj=8.2)
			
		}
	}
}











#creating PNG
tests<-nrow(expressionset) * ncol(pData(expressionset))
numberOfRows <- ceiling(tests / 3)
if(numberOfRows > 7){
	numberOfRows<-7	
}

png(filename = "../Images/genes_vs_clinical.png", width = 1440, height = 580*numberOfRows)
par(omi=c(0.2,0.2,0.2,0.2),mai=c(3,1,1,1))
layout(matrix(1:(3*numberOfRows),ncol=3,byrow=TRUE) )
plot_gene_vs_clinical(expressionset,probesets,label,annotation,horizontalScatterSpacing,doStatistics,verbose,cexFactoring=1.7)

if(numberOfRows > 7){
	frame()
	mtext("This is only a subset of available plots\nsee pdf for full set",cex=2)
}

garbage<-dev.off()







#creating PDF
pdf("../PDF/genes_vs_clinical.pdf", width = 3*5, height = 5*5)
layout(matrix(1:15,ncol=3,byrow=TRUE) )
plot_gene_vs_clinical(expressionset,probesets,label,annotation,horizontalScatterSpacing,doStatistics,verbose)
garbage<-dev.off()







#CREATING TXT
pdata<-cbind(data.frame(genesymbol=rep("",ncol(pData(expressionset))),genename=rep("",ncol(pData(expressionset)))), t(pData(expressionset)),stringsAsFactors=FALSE)
expression<-cbind(annotation[featureNames(expressionset),],as.data.frame(exprs(expressionset)))
rawData<-rbind(
		pdata,
		expression
)
colnames(rawData)[1:2]<-c("","")
write.table(rawData, file="../TXT/raw_data.xls", sep="\t",col.names=NA)
